%
% kernel.pro
%

:- module(kernel,
      [
        kernel/4,                        % +Kernel,+X,+Y,-Real
        kernel_modifier/5,               % +Modifier, +Kernel,+X,+Y,-Real

        % functions on lists
        list_to_multiset/2,              % +List,-Multiset
        cross_product/4,                 % +Predicate,+Xs,+Ys,-PXYs
        cross_product_write/3,           % +Predicate,+Xs,+Ys
        symmetric_cross_product_write/2, % +Predicate,+Xs
        common_stem/3,                   % +L1,+L2,-Len

        % functions on vectors
        inner_product/3,                 % +As,+Bs,-N
        multiset_inner_product/3,        % +As,+Bs,-N
        norm/2,                          % +As,-N
        diff/3,                          % +As,+Bs,-Cs

        % distances from kernels
        distance/4                       % +Kernel,+X1,+X2,-D

      ]).


% consult configuration file
:- (getenv('KERNEL_CONFIG',C) ; C='../config'), [C].


:- multifile
        default_kernel/2,                % +Type,-Kernel
        default_modifier/2.              % +Type,-Modifier

:- discontiguous
        kernel/4,
        kernel_parameter/3.

:- dynamic
        cross_product_cache/3.           % X,Y,PXY




%--------------------------------------------------------------------------
% TRACE KERNEL
%--------------------------------------------------------------------------

% comment or uncomment this kernel to trace calls to stdout
% NB. *must* be first kernel predicate encountered
%     or will miss values and leave unwanted choice points
/*
kernel(Kernel, X1::T1, X2::T2, 0) :-
    format('K_~w(\n\t~w :: ~w,\n\t~w :: ~w )\n', [Kernel, X1,T1, X2,T2]),
    fail.
*/



%--------------------------------------------------------------------------
% FUNCTIONS ON LISTS
%--------------------------------------------------------------------------

% NOTE: lists are represented here as native Prolog lists


%   list_to_multiset(+List, -Multiset)
%
%   converts prolog list to a multiset prolog list of value-multiplicity pairs
%   e.g. [b,a,b,c,c,b,c,c] -> [a-1,b-3,c-4]
%
list_to_multiset([], []).
list_to_multiset(List, Multiset) :-
    List \== [],
    msort(List, [H|T]),
    list_to_multiset(T, H-1, [], Multiset).

list_to_multiset([], H-Acc, MAcc, Multiset) :-
    append(MAcc, [H-Acc], Multiset).

list_to_multiset([H|T], H-Acc, MAcc, Multiset) :-
    Acc2 is Acc + 1,
    list_to_multiset(T, H-Acc2, MAcc, Multiset).

list_to_multiset([NewH|T], H-Acc, MAcc, Multiset) :-
    NewH \== H,
    append(MAcc, [H-Acc], MAcc2),
    list_to_multiset(T, NewH-1, MAcc2, Multiset).


%   cross_product(+P,+Xs,+Ys,-PXYs)
%
%   PXYs is the cartesian product of applying predicate P(+X,+Y,-PXY)
%
cross_product(P,Xs,Ys,PXYs) :-
    cross_product1(Xs,Ys,PXYs,P).

cross_product1([],_,[],_).
cross_product1([X|Xs],Ys,PXYs0,P) :-
    cross_product2(Ys,X,P,PXYs0,PXYs1),
    % writeln(PXYs0),
    cross_product1(Xs,Ys,PXYs1,P).

cross_product2([],_,_) --> [].
cross_product2([Y|Ys],X,P) -->
    { call(P,X,Y,PXY) },
    [PXY],
    cross_product2(Ys,X,P).


%   cross_product_write(+P,+Xs,+Ys)
%
%   writes the cartesian product of applying predicate P(+X,+Y,-PXY)
%
cross_product_write(P,Xs,Ys) :-
    cross_product_write1(Xs,Ys,P).

cross_product_write1([],_,_).
cross_product_write1([X|Xs],Ys,P) :-
    cross_product_write2(Ys,X,P),
    writeln(''),
    cross_product_write1(Xs,Ys,P).

cross_product_write2([],_,_).
cross_product_write2([Y|Ys],X,P) :-
    call(P,X,Y,PXY),
    write(PXY),
    (Ys \== [] -> write(',');true),
    cross_product_write2(Ys,X,P).


%   symmetric_cross_product_write(+P,+Xs)
%
%   writes the cartesian self product of applying predicate P(+X,+Y,-PXY)
%   optimised to only calculate the upper triangle using symmetry of P
%   Outputs squareform vector (inclusive of diagonal, which may be non-zero)
%
symmetric_cross_product_write(P,Xs) :-
    symmetric_cross_product_write1(Xs,P).

symmetric_cross_product_write1([],_).
symmetric_cross_product_write1([X|Xs],P) :-
    symmetric_cross_product_write2([X|Xs],X,P),
    symmetric_cross_product_write1(Xs,P).

symmetric_cross_product_write2([],_,_).
symmetric_cross_product_write2([Y|Ys],X,P) :-
    call(P,X,Y,PXY),
    write(PXY),
    write(' '),
    symmetric_cross_product_write2(Ys,X,P).


%   common_stem(+L1, +L2, -Len)
%
%   true when Len is the length of the common left sublist
%
common_stem(L1, L2, Len) :-
    common_stem(L1, L2, 0, Len).

common_stem([], _, Acc, Acc) :- !.
common_stem(_, [], Acc, Acc) :- !.
common_stem([H1|_], [H2|_], Acc, Acc) :-
    H1 \== H2.
common_stem([H1|T1], [H2|T2], Acc, Len) :-
    H1 == H2,
    NewAcc is Acc + 1,
    common_stem(T1, T2, NewAcc, Len).



%--------------------------------------------------------------------------
% FUNCTIONS ON VECTORS
%--------------------------------------------------------------------------

% NOTE: vectors are represented here as native Prolog lists of reals


%   inner_product(+A,+B,-N)
%
%   N is the inner product of vectors A and B
%
inner_product(As,Bs,N) :-
    length(As,L),                %mmm? checking length first may slow this down
    length(Bs,L),
    inner_product(As,Bs,0,N).
inner_product([],[],N,N).
inner_product([A|As],[B|Bs],N,Z) :-
    N1 is N + (A * B),
    inner_product(As,Bs,N1,Z).


%   multiset_inner_product(+A,+B,-N)
%
%   N is the inner product of vectors A and B
%
multiset_inner_product(A, B, N) :-
    multiset_inner_product(A, B, 0, N).

multiset_inner_product([], _, N, N).
multiset_inner_product([A|As], B, N, Z) :-
    multiset_inner_product2(A, B, N, N1),
    multiset_inner_product(As, B, N1, Z).

multiset_inner_product2(_, [], N, N).
multiset_inner_product2(Av-An, [Av-Bn|Bs], N, Z) :-
    N1 is N + (An * Bn),
    multiset_inner_product2(Av-An, Bs, N1, Z).
multiset_inner_product2(Av-An, [Bv-_|Bs], N, Z) :-
    Av \== Bv,
    multiset_inner_product2(Av-An, Bs, N, Z).


%   norm(+A,-N)
%
%   N is the norm of vector A
%
norm(A,N) :-
    inner_product(A,A,N1),
    N is sqrt(N1).


%   diff(+A,+B,-C)
%
%   C is the vector difference of vectors A and B
%
diff([],[],[]).
diff([A|As],[B|Bs],[C|Cs]) :-
    C is A - B,
    diff(As,Bs,Cs).

/*
% alternative version of diff/3 for adaptation: (probably slower though)
pdiff(X,Y,Z) :- Z is X - Y.
diff2(A,B,C) :- maplist(pdiff,A,B,C).
*/



%--------------------------------------------------------------------------
% DISTANCES FROM KERNELS
%--------------------------------------------------------------------------


%   distance(+Kernel, +X1, +X2, -D)
%
%   calculates euclidean distance D from any valid Kernel
%
distance(K, X1,X2, D) :-
    kernel(K, X1,X1, K11),
    kernel(K, X2,X2, K22),
    kernel(K, X1,X2, K12),
    D is sqrt(K11 - 2*K12 + K22).
%    D is sqrt(abs(K11 - 2*K12 + K22)).
    %FIXME: not sure why need abs here as -ve would break triangle inequality!!



%--------------------------------------------------------------------------
% KERNELS ON ATOMIC TYPES
%--------------------------------------------------------------------------


%   kernel(matching, +X1::Type, +X2::Type, -K)
%
%   K is defined for discrete data as 1 if X1=X2 else 0
%
kernel(matching, X1::Type, X1::Type, 1).
kernel(matching, X1::Type, X2::Type, 0) :-
    X1 \== X2.


%   kernel(product, +X1::{nat|int|real}, +X2::{nat|int|real}, -K)
%
%   K is defined for real data as |X1*X2|
%
kernel(product, X1::nat, X2::nat, S) :-
    S is X1 * X2.
kernel(product, X1::int, X2::int, S) :-
    S is abs(X1 * X2).
kernel(product, X1::real, X2::real, S) :-
    S is abs(X1 * X2).


%   kernel(difference, +X1::{nat|int|real}, +X2::{nat|int|real}, -K)
%
%   K is defined for real data as |X1-X2|
%
kernel(difference, X1::nat, X2::nat, S) :-
    S is abs(X1 - X2).
kernel(difference, X1::int, X2::int, S) :-
    S is abs(X1 - X2).
kernel(difference, X1::real, X2::real, S) :-
    S is abs(X1 - X2).



%--------------------------------------------------------------------------
% KERNELS ON VECTORS OF REALS
%--------------------------------------------------------------------------


%   kernel(linear, +X1::vector, +X2::vector, -K)
%
%   K is calculated as the difference between
%
kernel(linear, X1::vector, X2::vector, K) :-
    inner_product(X1,X2,K).


%   kernel(normalised_linear, +X1::vector, +X2::vector, -K)
%
%   K is calculated as the difference between
%
kernel(normalised_linear, X1::vector, X2::vector, K) :-
    inner_product(X1,X2,IP),
    norm(X1,N1),
    norm(X2,N2),
    N12 is N1 * N2,
    (N12 > 0    %FIXME: may need a very small no. instead of 0 as the test
    ->  K is IP / sqrt(N12)
    ;   K is 0
    ).


%   kernel(polynomial(L,P), +X1::vector, +X2::vector, -K)
%
%   K is calculated as (<X1,X2>+L)^P
%   where L is positive real and P is a positive natural
%
kernel(polynomial, X1::vector, X2::vector, K) :-
    % default values for L and R
    kernel(polynomial(1.0, 0.6), X1::vector, X2::vector, K).

kernel(polynomial(L,P), X1::vector, X2::vector, K) :-
    inner_product(X1,X2,IP),
    K is (IP + L) ** P.


%   kernel(gaussian_rbf(G), +X1::vector, +X2::vector, -K)
%
%   K is calculated as e^(-G * ||X1 - X2||^2)
%   where G is real
%   Note: this assumes euclidean space although the method used for
%         kernel_modifiers(gaussian(G),..) could be used instead here
%

kernel(gaussian_rbf, X1::vector, X2::vector, K) :-
    % default value for G
    kernel(gaussian_rbf(1), X1::vector, X2::vector, K).

kernel(gaussian_rbf(G), X1::vector, X2::vector, K) :-
    diff(X1,X2,D),
    norm(D,N),
    K is exp(-G * (N*N)).



%--------------------------------------------------------------------------
% KERNEL MODIFIERS
%--------------------------------------------------------------------------


%   kernel_modifier(default, +Kernel, +X1, +X2, -K)
%
kernel_modifier(default, Kernel, X1, X2, K) :-
    kernel(Kernel, X1, X2, K).


%   kernel_modifier(gaussian(+G), +Kernel, +X1, +X2, -K)
%
%   K is calculated as e^( -G * (K(X1,X1) -2K(X1,X2) + K(X2,X2)) )
%   where G is positive real and G > 0
%
kernel_modifier(gaussian(G), Kernel, X1, X2, K) :-
    kernel(Kernel, X1, X1, KX11),
    kernel(Kernel, X1, X2, KX12),
    kernel(Kernel, X2, X2, KX22),
    K is exp(-G * (KX11 -2*KX12 + KX22)).


%   kernel_modifier(normalised, +Kernel, +X1, +X2, -K)
%
%   K is calculated as K(X1,X2) / (K(X1,X1).K(X2,X2))^0.5
%
kernel_modifier(normalised, Kernel, X1, X2, K) :-
    kernel(Kernel, X1, X1, KX11),
    kernel(Kernel, X1, X2, KX12),
    kernel(Kernel, X2, X2, KX22),
    N12 is KX11 * KX22,
    (N12 > 0    %FIXME: may need a very small no. instead of 0 as the test
    ->  K is KX12 / sqrt(N12)
    ;   K is 0
    ).


%   kernel_modifier(polynomial(+P, +L), +Kernel, +X1, +X2, -K)
%
%   K is calculated as (K(X1,X2) + L)^P
%   where L is a positive real and P is a positive natural
%
kernel_modifier(polynomial(P, L), Kernel, X1, X2, K) :-
    kernel(Kernel, X1, X2, K0),
    K is (K0 + L)^P.



%--------------------------------------------------------------------------
% KERNELS FOR SIMPLE STRUCTURED DATA
%--------------------------------------------------------------------------

% "simple" here means that members of structures are homogeneous and atoms


%   kernel(list, +X1::list, +X2::list, -K)
%
%   K calculated by applying the matching kernel to elements on the lists
%   (which corresponds to computing p+q+r, where p is the length of the
%    shorter list, q is the number of consecutive matching elements at the
%    start of both lists, and r is 1 if the lists are the same length
%    otherwise r is 0)
%
%   X1,X2 are represented as Prolog lists  e.g. [a,b,c,b,a]::list
%
kernel(list, X1::list, X2::list, K) :-
    length(X1, LengthX1),
    length(X2, LengthX2),
    MinLength is min(LengthX1, LengthX2),
    common_stem(X1, X2, LeftMatch),
    (LengthX1 =:= LengthX2
     -> SameLength = 1
     ;  SameLength = 0
    ),
    K is MinLength + LeftMatch + SameLength.


%   kernel(set, +X1::set, +X2::set, -K)
%
%   K calculated by applying the matching kernel to elements on the sets
%   (which corresponds to computing the cardinality of the intersection or
%    to the inner product of the bit vectors representing the two sets)
%
%   X1,X2 are represented as Prolog lists  e.g. [a,b,c]::set
%
kernel(set, X1::set, X2::set, K) :-
    intersection(X1, X2, Intersection),
    length(Intersection, K).


%   kernel(multiset, +X1::multiset, +X2::multiset, -K)
%
%   K calculated by applying the product kernel to elements on the multisets
%   (which corresponds to computing the cardinality of the intersection or
%    to the inner product of the multiplicity vectors representing the two
%    multisets)
%
%   X1,X2 are represented as Prolog lists  e.g. [a,a,b,c,c,c]::multiset
%
kernel(multiset, X1::multiset, X2::multiset, K) :-
    list_to_multiset(X1, M1),
    list_to_multiset(X2, M2),
    multiset_inner_product(M1, M2, K).


%   kernel(multiset2, +X1::multiset2, +X2::multiset2, -K)
%
%   K calculated by applying the product kernel to elements on the multisets
%   (which corresponds to computing the cardinality of the intersection or
%    to the inner product of the multiplicity vectors representing the two
%    multisets)
%
%   X1,X2 are represented as Prolog lists of value-multiplicity pairs
%    e.g. [a-2,b-1,c-3]::multiset
%
kernel(multiset2, X1::multiset2, X2::multiset2, K) :-
    multiset_inner_product(X1, X2, K).


%   kernel(p_spectrum_binary(P), +X1::string, +X2::string, -K)
%
%   P is a positive natural
%
%   repetition insensitive version of p_spectrum
%   (ie. ignores no. times a match occurs in the strings)
%
%   X1,X2 are represented as Prolog atoms (not strings)
%    e.g. 'the cat'::string
%
%   (Note: this is my own accidental invention which arose
%    from an initial misreading of the spec in Nello's book)
%
kernel(p_spectrum_binary(P), X1::string, X2::string, K) :-
    setof(S, B^A^sub_atom(X1, B, P, A, S), S1),
    setof(S, B^A^sub_atom(X2, B, P, A, S), S2),
    intersection(S1, S2, Intersection),
    length(Intersection, K),
    !.
kernel(p_spectrum_binary(_), _::string, _::string, 0).
    % handle case when length of either X1 or X2 is less than P


%   kernel(p_spectrum(P), +X1::string, +X2::string, -K)
%
%   P is a positive natural
%
%   X1,X2 are represented as Prolog atoms (not strings)
%    e.g. 'the cat'::string
%
kernel(p_spectrum(P), X1::string, X2::string, K) :-
    bagof(S, B^A^sub_atom(X1, B, P, A, S), S1),
    bagof(S, B^A^sub_atom(X2, B, P, A, S), S2),
    intersection(S1, S2, Intersection),
    % WARNING: works because SWI intersection/3 allows duplicates
    %          in both input and output 'sets'
    length(Intersection, K),
    !.
kernel(p_spectrum(_), _::string, _::string, 0).
    % handle case when length of either X1 or X2 is less than P



%--------------------------------------------------------------------------
% KERNEL ON BASIC TERMS
%--------------------------------------------------------------------------


%   kernel(basic_terms, +X1::Type, +X2::Type, -R)
%
%   k : B x B --> R where B is a basic term and R is real
%
kernel(basic_terms, (X1::T)::Talias, (X2::T)::Talias, R) :-
    !,
    % Talias is a type alias for T so only affects modifier
    %%%FIXME: ignoring possibility of specific modifier on Talias for now
    kernel(basic_terms, X1::T, X2::T, R).

kernel(basic_terms, X1::T, X2::T, R) :-
    X1 =.. [C|L1],
    X2 =.. [D|L2],
    select_modifier(T, M),
    kernel_basic_terms(C,L1, D,L2, T, M, R).



%   kernel_basic_terms(
%       +TermType1, +TermBodyAsList1,
%       +TermType2, +TermBodyAsList2,
%       +KernelModifier,
%       -R
%   )
%
%   There are four cases in the definition of the kernel for basic terms:
%       case 1 - constants
%       case 2 - abstractions
%       case 3 - tuples
%       case 4 - non comparable
%   These four cases are dealt with below in a procedurally optimal order
%   rather than the order in which they appear in the kernel definition.
%
%   case 2 - abstractions: s,t in B_a where a = b --> g
%   k(s,t) = sum_(u in supp(s),v in supp(t)) k(V(s u),V(t v)).k(u,v)
%
kernel_basic_terms(abstraction,L1, abstraction,L2, _T, M, R) :- !,
    kernel_modifier(M, sum_kbt_over_cross_product, L1, L2, R),
    !.   % runs out of stack without this cut

%   case 3 - tuples: s,t in B_a where a = a_1 x ... x a_n
%   k(s,t) = sum_(i=1..n) k(s_i,t_i)
%            where s is (s_1,...,s_n) and t is (t_1,...,t_n)
%
kernel_basic_terms(tuple,L1, tuple,L2, _T, M, R) :- !,
    kernel_modifier(M, sum_kbt_over_lists, L1, L2, R).

%   case 1.1 - constants: C != D and s,t in B_a where a = T a_1...a_k
%   k(s,t) = K_T(C,D)
%
kernel_basic_terms(C,_L1, D,_L2, T, M, R) :-
    C \== D,
    !,
    select_kernel(T, KT),
    kernel_modifier(M, KT, C::T, D::T, R).

%   case 1.2 - constants: C = D and s,t in B_a where a = T a_1...a_k
%   k(s,t) = K_T(C,C) + sum_(i=1..n) k(s_i,t_i)
%            where s is C s_1...s_n and D is t_1...t_m
%
kernel_basic_terms(C,L1, D,L2, T, M, R) :-
    C == D,
    !,
    % RC = K_T(C,C)
    select_kernel(T, KT),
    kernel_modifier(M, KT, C::T, D::T, RC),
    % RP = sum_(i=1..n) k(s_i,t_i)
    kernel_modifier(M, sum_kbt_over_lists, L1, L2, RP),
    % R = K_T(C,C) + sum_(i=1..n) k(s_i,t_i)
    R is RC + RP.

%   case 4 - not comparable: a is not in universal set of
%            closed terms such that s,t in B_a
%   k(s,t) = 0
%
kernel_basic_terms(_,_, _,_, _, _, 0).



%   kernel(sum_kbt_over_lists, +Xs,+Ys, -R)
%
%   R is the sum of applying kbt to the lists Xs and Ys
%   where the elements of the lists are value::type
%
kernel(sum_kbt_over_lists, Xs,Ys, R) :-
    sum_kbt_over_lists(Xs,Ys,0,R).

sum_kbt_over_lists([],_,R,R) :- !.
sum_kbt_over_lists(_,[],R,R) :- !.
sum_kbt_over_lists([X|Xs],[Y|Ys],Acc,R) :-
    kernel(basic_terms, X,Y, R1),
    Acc2 is Acc + R1,
    sum_kbt_over_lists(Xs,Ys,Acc2,R).


%   kernel(sum_kbt_over_cross_product, +Xs,+Ys,-R)
%
%   R is the sum of kbt(V(s u),V(t v)) * kbt(u,v) to the cartesian product
%   Xs x Ys, where Xs is a lists whose elements are of the infix form
%   (s::type1)-(u::type2) and Ys is a lists whose elements are of the form
%   (t::type1)-(v::type2). The lists Xs and Ys represent sets here.
%
kernel(sum_kbt_over_cross_product, Xs,Ys, R) :-
    sum_kbt_over_cross_product(Xs,Ys,0,R).

sum_kbt_over_cross_product([],_,R,R).
sum_kbt_over_cross_product([X|Xs],Ys,Acc,R) :-
    sum_kbt_over_cross_product2(X,Ys, 0,R2),
    Acc2 is Acc + R2,
    sum_kbt_over_cross_product(Xs,Ys,Acc2,R).

sum_kbt_over_cross_product2(_,[],R,R).
sum_kbt_over_cross_product2(X-NX,[Y-NY|Ys],Acc,R) :-
    kernel(basic_terms, X,Y, R2),
    kernel(basic_terms, NX,NY, R3),
    Acc2 is Acc + R2 * R3,
    sum_kbt_over_cross_product2(X-NX,Ys,Acc2,R).



%--------------------------------------------------------------------------
% KERNEL SELECTION AND DEFAULTS
%--------------------------------------------------------------------------


%   select_kernel(+Type, -Kernel)
%
%   If a multifile default_kernel/2 predicate is defined for Type then succeed.
%   Otherwise, succeed with the locally defined default for type.
%   Only succeeds once.
%
select_kernel(Type, Kernel) :-
    default_kernel(Type, Kernel),
    !.
select_kernel(int,    product) :- !.
select_kernel(nat,    product) :- !.
select_kernel(real,   product) :- !.
select_kernel(_,      matching).


%   default_kernel(+Type, -Kernel).
%
%   Is true when Kernel is the default kernel for Type
%
%   Multifile predicate tested by select_kernel/2 to allow
%   application-specific declaration of kernel bindings
%
    %NB. no definition is required in the current module


%   select_modifier(?Type, ?Modifier)
%
%   If a multifile default_modifier/2 predicate is defined for Type then succeed.
%   Otherwise, succeed with the default modifier.
%   Only succeeds once.
%
select_modifier(Type, Modifier) :-
    default_modifier(Type, Modifier),
    !.
select_modifier(_, default).


%   default_modifier(?Type, ?Modifier).
%
%   Is true when Modifier is the default kernel for Type
%
%   Multifile predicate tested by select_modifier/2 to allow
%   application-specific declaration of kernel modifier bindings
%
    %NB. no definition is required in the current module



